function cook_membrane
clc
%close all
addpath('../');
addpath('../../');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% create mesh/read mesh
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% projDir = 'home/hbui/kratos_janosch';
% projName = 'cook-membrane.gid';
%fn = [projDir '/' projName '/' 'mesh4.msh'];
fn = '/home/hbui/kratos_janosch/cook-membrane.gid/mesh4.msh';
mg = fem2.gid.GidMeshGenerator(fn);
mesh = mg.getMesh(2);
M = fem2.Model(mesh);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% conditions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
meshUtil = fem2.MeshUtilities(mesh);
c1 = fem2.Constraint([0 0]); %fix x, fix y
leftedge = meshUtil.seekNodesOnSurface(1,0,0);
rightedge = meshUtil.seekNodesOnSurface(1,0,-48);
tip = meshUtil.seekNode([48 60]);
mesh.setConstraint(c1,leftedge);
F = 100;
for i=1:rightedge.size()-1
    ldx = math2.ConstantFunctionRnToR(0);
    ldy = math2.ConstantFunctionRnToR(F);
    n1 = rightedge.get(i-1);
    n2 = rightedge.get(i);
    e = mesh.addEdge(fem2.Edge([n1 n2]));
    el = fem2.EdgeLoad(e,[ldx ldy]);
    M.addLoad(el);
end
% for i=1:rightedge.size()
%     n = rightedge.get(i-1);
%     f = fem2.Force(n,[0 F]);
%     M.addLoad(f);
% end
mesh.setNumberOfFields(2);
plotMesh(mesh,{'normal','conditions','id'});
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% element declaration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nf = mesh.countFaces();
nu = 0.4999;
E = 250;
C = C_plane_strain(E,nu);
for i=1:nf
    e = fem2.MatlabStructuralElement(mesh.getFace(i-1));
%     e.setIntegrationOrder(2);
    M.addElement(e);
    es = createElementContainer();
    es.constitutiveMatrix = C;
    elements(i) = es;
end
elements = Initialize(elements,M);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
postProcessInfo.infoList = {'resistant force',...
    'stress at Gauss points'};
% postProcessInfo.infoParam = {[1 4],[1 2]};
postProcessInfo.infoParam = {[1 4],[1]};

lambda = 1;
[elements,output] = solving_strategy_static_2(elements,M,...
    lambda,50,1e-6,postProcessInfo);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% post process
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% res_node_1(:,:) = -output.resistant_force(:,1,:)
% delta = lambda*du;
% plot(delta,res_node_1);

plotDeformedMesh(mesh,1);

disp(['number of dofs: ' num2str(mesh.countDofs())]);
disp(['tip displacement: ' num2str(tip.getUHatLocal()')]);

rmpath('../');
rmpath('../../');
end

